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Abstract 

Many features of the sequence of action potentials produced by repeated stimulation of a cardiac patch can be modeled 
by a ID mapping, but not the full behavior observed in the restitution portrait: in particular, not (i) distinct slopes 
for dynamic and S1-S2 restitution (rate dependence) and not (ii) long transients in the approach to steady state 
(accomodation). To address these shortcomings, ad hoc 2D mappings, where the second variable is a "memory" 
variable, have been proposed; it seems that these models exhibit some, but not all, of the relevant behavior. In 
this paper we introduce a new 2D mapping and determine a set of parameters for it that gives a rather accurate 
description of the full restitution portrait found for one animal. The fitting procedure is straightforward and can 
easily be applied to obtain a mathematical model for data from other experiments, including experiments on different 
species. Uniqueness of the parameter choice is also discussed. 

The changes in the mapping, compared to previous models, result from requiring that the mapping can be derived 
as an asymptotic limit of a simple ionic model. Among other benefits, one can interpret the parameters in the mapping 
in terms of the ionic model. The ionic model is an extension of a two-current model that adds a third dependent 
variable, a generalized concentration. The simplicity of the ionic model and the physiological basis for the mapping 
contribute to the usefulness of these ideas for describing restitution data in a variety of contexts. 



1 Introduction 

Sudden cardiac death kills half of all people who die of heart disease, the Number One cause of death 
in the United States, accounting for more than 400,000 deaths each year As part of a larger effort 
to understand and prevent arrythmias that lead to death, much work has been devoted to studying the 
resitution properties of cardiac tissue under repeated stimulation, especially rapid stimulation EH El E] • 
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The restitution portrait, as is obtained in the recently proposed perturbed downsweep protocol |Sj , combines 
data about both dynamic restitution and S1-S2 restitution, as well as information about the transients in 
the approach to steady state. (These terms are defined in the Appendix.) Because the restitution portrait 
includes such a wealth of information, experimental data from this protocol provide a demanding test of 
mathematical models for restitution. Specifically, restitution portraits obtained from small pieces of frog 
myocardium jS] show two noteworthy behaviors: 

1. Under rapid pacing but still exhibiting a 1:1 response, the dynamic restitution curve is typically quite 
steep, having slope 3 or even greater, while the S1-S2 restitution curve is rather shallow, having a slope 
on the order of 0.5. 

2. The approach to steady state is very slow, with a time constant on the order of 20-30 sec. The 
approach, which is monotonic, occurs along a smooth curve that is distinct from either restitution 
curve. This behavior is observed over the entire range of pacing intervals leading to a 1:1 response. 

These behaviors illustrate phenomena in cardiac electrophysiology variously known as rate dependence, mem- 
ory, and accomodation. 

Consider a typical experiment of producing a sequence of action potentials by the application of periodic 
stimuli to a small patch of cardiac tissue 1 . Let B (mnemonic for basic cycle length) be the period of the 
stimulus application, and let A n be the duration of the n th action potential. Nolasco and Dahlen 7 proposed 
to describe the sequence A n approximately by a recursion relation 

A n+1 = G(B - A n ) (1) 

for an appropriate function G. Of course B — A n is the diastolic interval D n following the n th action potential. 

Although the introduction of ID mapping models was a decisive advance in our understanding of resti- 
tution, such models cannot exhibit Behaviors 1 or 2 above; indeed, for a mapping of the form |JI|I. dynamic 
restitution, S1-S2 restitution, and the approach to steady state all fall on the same curve. In an effort to 
model these behaviors, various authors [SJIH] have introduced ad hoc 2D mapping models, one variable being 
A n and the other a memory variable. In unpublished analysis we have shown that parameters may be chosen 
in these models to obtain some of the behavior in [Jy, but we believe that no single set of parameter values 
can produce the complete range of behavior. 

In this paper we introduce a 2D map that 

(i) fits all the dynamic restitution data, all the S1-S2 restitution data, and most of the transient data from 
the restitution portrait 6 that involves a 1:1 response (in particular, exhibiting both Behavior 1 and 
Behavior 2) 

and moreover 

(ii) can be derived with asymptotic analysis from a simple ionic model. 

Regarding (i), determining the best parameter fit is a straightforward process, and it may be applied to 
model other restitution experiments with different animals. Regarding (ii), in this ionic model, memory is 
associated with the accumulation over time of ions in the cell. The ionic model suggests some subtle changes 
in the mapping compared to previous models of this type |SJ |U] , and these changes led to a good fit with 
experiment. 

1 Of course the heart is a spatially extended structure. In performing experiments with a small piece of tissue, one hopes 
that propagation effects can be eliminated, in order to isolate a simpler problem. In any event, in this paper we consider models 
without propagation. 
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2 The mapping model and its ionic basis 



Like 8,9, our mapping contains two variables — one is A n as above; the other, which will be denoted by c„, 
specifies the ionic concentration in the underlying ODE at the start of the n th action potential. We suppose 
that, under periodic stimulation (say with period B), these two variables evolve according to the iteration 

A n+1 = G(D n ) + $(c„ +1 ) (2) 

c n+1 = (c n + e )e- B ' T *™* (3) 

where 

G(D) = A ma * + r fclosc In { 1 - ae- D ' T °^ } (4) 

and 

*(c)=r_ln{^}. (5) 

In J2HU), A max , T fclOBC , t sc1oso , T opon , and r pump are parameters with the dimensions of time, and e,a, and (3 are 
dimcnsionless parameters. As we discuss below, this mapping was derived as an asymptotic limit of an ionic 
model, and the constants in (0] [SJ refer to constants in the ionic model. Thus, by relating the mapping 
to the ionic model, we will be able to understand the significance of these constants, as well as explain the 
mnemonics of their subscripts. This physiological connection enhances the utility of the mapping in modeling 
various experiments. 

The mapping 10 EI) must be supplemented with the restriction that no action potential is produced unless 
the diastolic interval (DI) is long enough. Physically, if the DI is too short, the system will skip a beat and 
may jump to a 2:1 response. Mathematically, if the DI is less than r opon In a, then is not even defined. The 
precise lower bound for the DI, which depends on the concentration c, cannot be derived without analysis 
of the ionic model. These somewhat technical issues will be addressed in a future publication. 

The ionic model underlying 1013) is a modest extension of the two-current model |l()llllj in that: (a) the 
equations for the two variables in the two-current model, the voltage v and the gate h, are slightly more 
complicated and (b) a third variable, the concentration c, is added to the system. Specifically: 

(i) In the modified equation for v, 

dv 

— = -J in (v, h, c) - J out («), (6) 

the outward current is still given by J out (w) = w/r out , but the inward current is now the sum of concentration- 
independent and concentration-dependent parts 

Mv, h, c) = - A {<j> ci (v) + e- c ^ d (v)} . (7) 

The response of the model to stimulation is not very sensitive to the exact form of the functions 4> ci (v) and 
(f> ad (v), and we do not specify these functions here. 

(ii) In the modified equation for h, 



dh 
~dt 



(l-h)/T open if v < v„ it 

(8) 

-h/T olaae (v) ifu>i; crit , 
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the time constant for the closing of the gate depends on the voltage v. The precise form of T closo (u) is less 
important than the fact that two different time-scale parameters derive from the closing of the gate, say 

r fo i ose : fast closing for v near the top of its range, and 

Tsdose'- slower closing for smaller v (but still greater than v am ). 

(iii) The equation 2 for c may be written 

dc 

— = I{t) - c/r pump . (9) 

Thus, the concentration is determined by a balance between lit), the current which leads to the build-up of 
charge in the cell, and linear pumping, which removes charge from the cell. Again, the precise form of I(t) 
is not important, only two key properties 3 : 

• I(t) is nonzero only during the upstroke of an action potential (i.e., the initial rapid increase of v), and 

[ I(t)dt = e, (10) 
J J 

where J denotes the interval between two consecutive stimuli and e is the same-named parameter in Q- In 
words, lfTU|l asserts that a fixed charge e enters the cell during each action potential. 
Under the assumption that 

T in < T om < T tcloBC ,T 5cloao ,T opcn < T pump , (11) 

the mapping 13 may be extracted as an asymptotic limit of 10 00 EH' ■ 

The significance of the parameters r pump , r fcloao , r scloso , and r opon in J3[S1) is clear from the discussion of the 
ionic model. A max is the longest possible action potential since in 1 1151) . the terms 

r Mose In { 1 - ae- D 'T<»~ } , r sclose In j j 

are nonpositive and tend to zero under repeated slow pacing. As noted above, e specifies the amount of charge 
that enters the cell during each action potential. The parameter a is a dimensionless combination of several 
parameters in the ionic model. Similarly, [3 is the ratio of the strengths of the concentration-independent 
and concentration-dependent parts of J in (v,h,c). 

3 Fitting the model to experimental data 

The model (|2H3) contains eight parameters. One might attempt to determine all these parameters simul- 
taneously with a massive least-squares fit to the entire restitution portrait. However, there probably are 

2 Equation J^J appears to be nonautonomous. There is an autonomous form of the equation in which I(t) is expressed in 
terms of the voltage v(t). However, I(t) is a function of dv/dt rather than v itself, and we prefer the notation in 1^1 . 

3 The assumption II 1)1 was motivated by some numerical experiments with the LRd model 1121 . Specifically, we found that 
continued rapid stimulation, which leads to shorter APD's, also increases the concentration of sodium ions in the cell. If the 
model is modified by artificially fixing the Na+ concentration at its resting value, the APD's are shortened to a lesser degree 
1131 . Condition I1UI is designed to mimic the fast sodium current, which flows at the start of the action potential and only as 
long as needed to depolarize the cell. Although such considerations motivated our assumptions, we regard the good fit of the 
model with experimental data as stronger supporting evidence. 
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multiple local minima of the residual in such a high dimensional parameter space. Even if there were a 
unique minimum, there may be a large set in parameter space that gives a fit of essentially the same quality, 
and in high dimensions it is difficult to assess this prospect quantitatively. Moreover, such a massive fit 
would not yield any insight about the significance of the individual parameters. 

The restitution portrait provides an attractive alternative — we determine parameters sequentially, in 
small groups, from four different types of data, as follows: 

(f) Slope of the S1-S2 restitution curve: r fc i oao ,a 

(2) Slope of the dynamic restitution curve: r sclOB( . , er pump , [3 

(3) Overall height of the dynamic restitution curve: A max . 

(4) Transient time constant: r pump 

Note that T opcn is not on the above list. As we shall see below, it turns out that even this sequential approach 
leaves r opon quite uncertain; one obtains only a modest lower bound for this parameter. 



3.1 S1-S2 restitution data 

To determine the slope S12 experimentally, many stimuli (say N) are applied with a fixed basic cycle length 
(say B) until the system achieves a steady state, after which the (N + l) sf stimulus is applied following a 
perturbed cycle length B + A: then S%2 is defined by 

_ A N+1 - A N 

bl2 ~ n n ■ " 2 > 

Note that 

A n = A ss for n large and n < N, 
D„ = D ss for n large and n < N, 

while 

D N = D SS + A, so that — -Djv-i = A. 

If, under the above protocol, a sequence of action potentials is determined by the mapping then 
the associated concentrations also satisfy 

c n = c ss for n large and n < N: 

moreover, since by JTTJ, r pump is so large, 

CN+l ~ Cat. 

Thus, the model J5] |3J) predicts that 

S12 w G'(D SS ). (13) 
Differentiating Q, we see that the RHS of (fH5|) is given by 

G'(D SS ) = Tfc '- D 1 . (14) 

The squares in Figure I show S12 as obtained in one of the experiments of [o]. It is natural to try to 
choose the three parameters in ljI4ll . r fclOBC , r opcn , and a, by finding the best fit to this data. However, the 
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Figure 1: A fit of HU) to S 12 data, with r opon = 300 msec 
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Figure 2: The residual in fitting (|14Jl to S12 data, as a function of r open 

quality of the fit is not at all sensitive to r opon , provided this parameter is sufficiently large. This is illustrated 
in Figure 2, which shows the residual as a function of r opcn when the corresponding optimal values of r fclOBO 
and a are chosen. Somwhat arbitrarily, in the remainder of the fitting process we have chosen a minimal 
value for this parameter, 

t opc „ = 300 msec. (15) 

The corresponding optimal values for r fcloao and a are given in Table 1, and the fit of the Sf-S2 data from 
the experiment for this value of r opcn is given in Figure 1. Of course, because the parameters are coupled 
in various ways, the ambiguity in r open means that most of the parameters cannot be regarded as uniquely 
determined by the experimental data. 



G 



Parameter 


Value 


Units 


T fc l ose 


29.4 


msec 


a 


1.174 




7"s close 


570 


msec 


FT 

' pump 


821 


msec 


P 


0.386 




T 

pump 


30000 


msec 


A max 


889 


msec 



Table 1: Optimal parameter values, assuming T opon = 300 msec 

3.2 Slope of the dynamic restitution curve 

According to (J2J), the steady-state action potential duration is given by 

A ss = G(D ss ) + $(c ss ). (16) 

In the physiological range, B is larger than, but on the order of, a typical action potential. The order of 
magnitude of the action potential duration in the model is determined by the gate constants r fclOBO , r acloao , r opcn . 
But by H 1 1 p . the gate constants are much smaller than r pump . Thus, we conclude that, in the physiological 
range, 

B«r pump . (17) 

Given ifPTjl. it follows from © that 

Css « er pump /B. (18) 

Combining (|16l I18f) , we find that 

A ss = G{D SS ) + $(er pump /fl). (19) 

Since B = A ss + D ss , (|19|l defines ^4 SS implicitly as a function of D ss , and (|19|l may be differentiated to 
yield the slope of the dynamic restitution curve 

Sdyn = -jjy- = G'{D SS ) - — (Sdyn + 1), (20) 

where we have used the fact that dB/D ss = Sd yn + 1- Recalling (|13|l . computing d&/dc from and 
rearranging, we obtain 

Sdyn <Sl2 7" a closc£ T pump 1 /„-. 

Sdyn + l & 1 + p e ^JB ■ \ > 

The squares in Figure 3 show data for the combination of the slopes Sd yn and £12 on the LHS of (|21|l as 
obtained in one of the experiments 4 of [S]. The curve in the figure shows the best least-squares fit 5 , and the 
corresponding parameter values for r BCloao , er pump , /? are listed in Table 1. 

4 More accurately, data for Sjj„ were taken from the experiment, but values for S12 from our fit of 1141 were used. 

5 Of the various parameter fits in the paper, Figure 3 is perhaps the least satisfactory. Note that this weakness is far less 
clear in the fit of Figure 4, an integrated version of 1211 . This difference illustrates the point that data from the restitution 
portrait provide a demanding test of any mathematical model. 
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Figure 3: A fit of (|21(l to dynamic restitution slopes 
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Figure 4: A fit of (0 El to the dynamic restitution curve 



3.3 The two remaining parameters 

We choose A maX) which sets the overall height of the restitution curve, to fit the data as shown in Figure 4. 

Steady state data determines the product er pump but neither factor individually; for the individual factors, 
we turn to the transient data. In the downsweep of the restitution-portrait protocol [H], each time the basic 
cycle length is decreased, a time constant is extracted from the transient to the next steady state. These 
time constants are typically on the order of 20-30 sec. On the other hand, in the mapping model EJ, 
the time constant of such transients is very nearly equal to r pump for all B. Although the time constants 
in the experiments vary somewhat with B, we merely take r pump = 30,000 msec, as indicated in Table 1, 
which yields e = 0.0274. (For a possible improvement of this part of the model, see the "Nonlinear pumping" 
discussion in Section 4.3.) 
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4 Discussion 



4.1 Conclusions 

In this paper we have introduced an ionically based 2D mapping model for cardiac restitution. The underlying 
ionic model suggests some subtle changes from previous models (B1E]j an d having made these changes, we are 
able to choose parameters in the mapping to fit S1-S2 restitution data, dynamic restitution data, and most 
transient data, as obtained in [S] . The simplicity of the procedure means that mathematical descriptions of 
various experiments could be obtained with these techniques. 

We obtained parameters in J3] 01 by fitting different types of data derived from the restitution portrait. 
The complete restitution portrait from the experiment [S] is shown in Figure 5(a). As explained in [B] (also 
see the Appendix), such a restitution portrait presents in compact form the dynamic restitution data, the 
S1-S2 restitution data, and the transient data obtained in the perturbed downsweep protocol. In Figure 5(b) 
we show the restitution portrait from the map. The agreement is remarkably good, with one caveat: In the 
experiment, the basic cycle length was decremented in steps of approximately 100 msec. In computations 
with the mapping, at the two shortest BCL's, steps this large led to initial diastolic intervals less than the 
value T opon In a where Q becomes undefined. To circumvent this difficulty, it was necessary to decrease the 
BCL in smaller steps. Because of this issue, at the two shortest BCL's, the transient shown for the mapping 
is shorter than for the experiment. (Remark: In the experiment, data from late in the transients is not 
recorded and hence is not shown in Figure 5(a).) 

In a nonlinear least-squares fitting of data, one wonders whether other, radically different, parameter 
values could fit the data equally well. For example, as we saw, the quality of the fit was spectacularly 
insensitive to T opcn . Accepting the ambiguity with respect to r opon , we further addressed this issue as follows: 
In our fitting of lfH|l to S12, we evaluated the 2x2 Hessian matrix of the residual with respect to T fclosc 
and a; all diagonal entries of this matrix were large; one of its eigenvalues was also quite large, but the 
other was rather small. The situation for fitting (|21(l was similar, with one even smaller eigenvalue. From 
this we conclude that, with r opon given, changing any one other parameter would affect the quality of the fit 
adversely, but groups of parameters might be changed with little effect. The issue is still under investigation, 
in particular how it might affect the transients at shorter basic cycle length. 

4.2 Comparison with previous models 

The memory variable M n in [HI IH) evolves according to 



This equation differs only slightly from J3J), but this small change was crucial in obtaining a model that 
exhibited both Behavior 1 and Behavior 2, which were seen in experiments. Equation 10) results from our 
assumption l|10|) that a fixed charge e enters the cell during each action potential. The evolution i|22|) can 
also be obtained from an ionic model: specifically, a model in which charge accumulates during the action 
potential and is pumped out during the diastolic interval. 
In E], the analogue of G(D) in (0) is given by 



where K\,K%, and r are constants. The shape of i|23|) is qualitatively similar to that of @, with K\,K^ 1 
and t in 123|) playing approximately the same role as ^4 max , T fclOBO , and r opon , respectively. However, (0) 




(22) 



-D/t 



(23) 
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(a) 




Figure 5: The full restitution portrait, (a) Experiment (b) The mapping J2HHJ) 



contains the fourth parameter a, which is crucial in fitting the S1-S2 restitution data of UJ. Incidentally, 
the logarithmic dependence in is what leads to the fact that the memory effect in is additive, not 
multiplicative as in Since Q is based on an ionic model, we regard the additive form as more natural, 

but either form could be fitted to the data. 

In 14 . Elharrar and Surawicz parametrize the dependence of the steady-state action potential A ss on 
the period B with a function that has the asymptotic form as B — > oo 

A ss ~ A max - Const/B. (24) 

This algebraic decay is important since it matches the experimental observation that ^4 SS continues to increase 
significantly as B is raised far beyond physiological values; by contrast, the exponential decay in (|23[) would 
not lead to this behavior. Let us show that J21I2J) exhibits such algebraic decay for large B. Since G decays 
exponentially, it will contribute a constant term to l|24|) but no term at the level 0(1/ B). Recalling that 
(f 1 71) holds in the physiological range and referring to |@| [5J , we obtain (|24|) on substitution into (|16|l , where 
the coefficient of 1/B equals 

7"sclose£7"pump 

1 + fi ' 



4.3 Future work 

(a) Alternans 

In this paper we have considered only 1:1 responses of the heart. Because of the model's transparancy, 
it is not difficult to adjust parameters in J21I3|) to produce a 2:2 response, or alternans. However, this does 
not seem to lead to a good model for the experiments [H], for the following reason: 

A mapping leads to alternans when an eigenvalue, say A, of its Jacobian at the fixed point passes through 
-1. Because r pump is so large, the concentration c has little influence on the relevant eigenvalue; specifically 

A a -S 12 . (25) 

However, it is common in the experiments to find alternans even when S12 is on the order of 0.5. We hope 
to derive a better model for alternans by including calcium fluxes, as in |15|. 



(b ) Nonlinear pumping 

In our ionic model, © describes a linear pump that pumps at a very slow rate. For this model, the time 
constant for the transient leading to steady state is very nearly equal to r pump , no matter what the basic 
cycle length B is. However, experiments (e.g., Figure 4 of ^3]) show that this time constant decreases if B 
is decreased. Our model could be modified in several ways to better reproduce this behavior. Perhaps the 
simplest is to allow nonlinear pumping: i.e., to change JHJ to 

1-/(0-*- 

for some power p. This will be developed in a future publication. 



(c) Global S1-S2 restitution 

In S1-S2 experiments, after steady state has been attained with stimuli spaced at an interval SI, the 
APD following a perturbed stimulus interval S2 is measured. By a global S1-S2 restitution curve, we mean 
the graph of action potential duration as a function of S2, when S2 varies over its entire range (with SI 
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fixed). By contrast, the restitution portrait |S] studies only what might be called local S1-S2 restitution: 
i.e., only perturbed intervals 52 nearly equal to 51 are explored. Of course global restitution provides more 
data to test a model, and experiments are planned to test the predictions of the model (J^I EJ m this regard. 

One positive conclusion may already be reported. Elharrar and Surawicz |14) found that each global 
S1-S2 curve crosses the dynamic restitution curve twice, both for 52 = 51 and again for small DI. Observe 
from JIJ that G(D) tends to negative infinity as D tends to r opcn In a > from above. Thus, the S1-S2 and 
dynamic restitution curves produced by J5] UJ also intersect twice in this way. 

Support of the National Institutes of Health under grant 1R01-HL-72831 and the National Science Foun- 
dation under grants PHY-0243584 and DMS-9983320 is gratefully acknowledged. 
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Figure 6: Schematic action potentials 



Appendix A 

In this Appendix, we briefly summarize the notation used throughout the paper and the protocols used to 
generate the various restitution curves. 



Tissue response 

When a small piece of cardiac muscle is subjected to a sequence of brief electrical stimuli who strength exceeds 
a critical threshold, the myocytes respond by producing action potentials. Figure 6 shows a schematic of the 
transmembrane voltage measured from a single myocyte. Following a stimulus, first the voltage rises rapidly 
(indicating cell depolarization) , then it has a plateau region during which the cell cannot be reactivated (the 
refractory period), and finally the voltage returns to its resting value (the cell repolarizcs). This time 
course is known as an action potential, and the duration of the n th action potential is denoted by A n . The 
interval between the time when the cell repolarizes (following the n th action potential) and the time when 
it depolarizes again (due to the (n + l) st stimulus) is known as the diastolic interval and is denoted by D n . 

For a periodic train of stimuli delivered at a slow rate, the myocytes display a phase-locked 1:1 response, 
where each stimulus produces an identical action potential. At faster pacing rates, the 1:1 response pattern 
is sometimes replaced by a 2:2 phase-locked period-2 response pattern, known as alternans. Under other 
conditions, either the 1:1 or 2:2 response becomes unstable and is replaced by a 2:1 response pattern, where 
only every other stimulus elicits an action potential. 



Electrical restitution 

A fundamental characteristic of cardiac cells is the shortening of A n as the pacing rate increases, known 
as electrical restitution. A restitution curve quantifies the functional relation between D n and A n . As 
pointed out by Elharrar and Surawicz the restitution curve depends on the pacing protocol used in its 
measurement. Several protocols are currently in use today, as described in detail in and briefly below. 

Dynamic restitution curve protocol 
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Pace the tissue periodically with period B until the tissue reaches steady state. Record A n+ i = A* and 
D n = D* . Change B and repeat. The pairs of points (D* , A*) form the dynamic restitution curve. 

S1-S2 restitution curve protocol 

Pace the tissue periodically with period B (the "SI" stimuli in "S1-S2") until the tissue reaches steady state. 
Make a sudden change in the pacing period (the "S2" stimulus), record the diastolic interval D n preceding 
the S2 stimulus and the resulting action potential duration A n+ \. Repeat for various S2 intervals. The 
pairs of points (D n , A n+ i) form one of the S1-S2 restitution curves. There is a different S1-S2 restitution 
curve for each value of SI, a characteristic of rate-dependent restitution. 

Restitution portrait 

The restitution portrait is a combination of both protocols 6 described above, where, for a given B, only 
segments of the S1-S2 restitution curves in a neighborhood of the corresponding steady state (D* ,A*) are 
measured. The restitution portrait also shows the pairs (D n , A n+ {) in the transient leading to each steady 
state. 



In this paper we do not consider the constant-B protocol of 



14 



